Non-perturbative approach to the attractive Hubbard model 



S. Allen, and A.-M.S. Tremblay 1 
Departement de Physique and Centre de recherche sur les 
proprietes electroniques de materiaux avances. 
1 Institut canadien de recherches avancees 
Universite de Sherbrooke, Sherbrooke, Quebec, Canada J1K 2R1. 
(February 1, 2008) 

o 

A non-perturbative approach to the single-band attractive Hubbard model is presented in the 
■ general context of functional derivative approaches to many-body theories. As in previous work on 

the repulsive model, the first step is based on a local-field type ansatz, on enforcement of the Pauli 
principle and a number of crucial sum-rules. The Mermin- Wagner theorem in two dimensions is 
automatically satisfied. At this level, two-particle self-consistency has been achieved. In the second 
step of the approximation, an improved expression for the self-energy is obtained by using the 
, results of the first step in an exact expression for the self-energy where the high- and low-frequency 

behaviors appear separately. The result is a cooperon-like formula. The required vertex corrections 
are included in this self-energy expression, as required by the absence of a Migdal theorem for 
. this problem. Other approaches to the attractive Hubbard model are critically compared. Physical 

,1 ^ 1 consequences of the present approach and agreement with Monte Carlo simulations are demonstrated 

, in the accompanying paper (following this one). 

^ ■ PACS numbers: 71.10.Fd, 71.27.+a, 71.10.-w, 05.30.Fk. 
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I. INTRODUCTION 



In the mid to late 1950's, quantum-field theoretical methods that had been developed first in the context of quantum 
T—t \ electrodynamics began to have widespread applications in Condensed Matter Physics. [1,2] One can, roughly speaking, 
^ • distinguish two types of approaches. The diagrammatic methods of Feynman and the functional methods of the 
ON ' Schwinger school [3,4]. Both points of view on many-body theory are equivalent. In particular, perturbation theory 
can be formulated diagrammatically or with functional methods. The two approaches in fact complement each other. 
For example, in calculating response functions, subsets of diagrams are often summed to infinite order. But naive 
resummations will generally break gauge invariance or other exact symmetry properties unless consistency between 
self-energy and two-particle irreducible vertices is enforced following a technique whose most natural formulation 
| employs functional derivatives [4,5]. 
"jS ■ Diagrammatic methods have nevertheless become by far the most popular techniques for many-body [6] problems 
in Condensed Matter, but the quest for non-perturbative approaches leads in general outside the realm of diagrams. 
While the Hartree-Fock approximation has a diagrammatic interpretation, what seems to be the most accurate 
approach to the electron gas at metallic densities (local field approximation (LFA) [7]) does not have a simple dia- 
grammatic interpretation. In the case of the Hartree-Fock approximation, a variational principle guides the accuracy 
of the approximation. In the LFA, it is a self-consistency requirement at the two particle level that controls the 
q | accuracy. 

Non-perturbative approaches have been developed in particular for the Hubbard model [8] , perhaps the best know 
. £h ' model for strongly interacting electrons on a lattice. While the early Green function decoupling schemes have largely 
fallen in disfavor because of their ad hoc and uncontrolled nature, many successful non-perturbative approaches 
have been developed in one dimension. Recently, dynamical mean-field theory has become the method of choice [9] 
for higher dimensions. This method, however, is in fact based on an expansion about infinite dimension. In two 
dimensions in particular, the momentum dependence of the self-energy cannot be neglected and this method becomes 
less accurate. 

The purpose of this paper is to extend to the attractive Hubbard model a non-perturbative approach developed 
previously for the repulsive model [10-12]. This method was an extension of the LFA work of Singwi, Tosi, Land and 
Sjolander [7] for the electron gas and of Hedeyati and Vignale for the Hubbard model [13]. It went further than these 
work in imposing the Pauli principle, a number of exact sum rules, of conservation laws, and proposing a formula 
for the self-energy in the paramagnetic state that includes momentum and frequency dependence. The approach 
also includes an internal check on accuracy based on an exact relationship between one and two-particle properties. 
Although the method fails in strong coupling or very close to a critical point, it gives the most accurate results 
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when compared with Monte Carlo simulations in the weak to intermediate coupling regime [12]. The present paper 
generalizes that approach to the attractive Hubbard model. The accompanying paper [14] demonstrates accuracy by 
comparisons with Monte Carlo simulations, and discusses a problem of importance in the context of high-temperature 
[15] and organic [16] superconductors, namely the opening of a pairing-fluctuation induced pseudogap. 

The structure of this paper is as follows. In section II we introduce the general many-body formalism to obtain 
expressions for the self-energy and irreducible vertices in the functional derivative approach. Although in principle 
standard, the functional derivative approach in the particle-particle channel is not widely used. It allows us to 
establish a number of exact results. At this level, all the results could have been obtained with formal diagrammatic 
expansions in skeleton diagrams, but we refrain from doing this since the functional approach is far more economical 
in this context. The exact results that we present form the basis of the Two-Particle Self- Consistent (TPSC) approach 
and its extension (section III). These arc non-perturbative approximations that are not diagrammatic resummations. 
A discussion of the various exact results that our approach satisfies either exactly or self-consistently is then presented 
in section IV. Details of the derivation of the exact results, namely sum rules and high-frequency expansions, arc 
given in Appendix A. The last appendix compares our approach with various other approaches and explains the 
connection with the formalism for the repulsive Hubbard model. A summary of our main results is in the conclusion, 
section V. 



II. EXACT RELATIONSHIPS BETWEEN GREEN FUNCTION, SELF-ENERGY AND VERTICES 

In this section, we derive a number of exact results using the functional derivative formalism. It is on the basis of 
these results, and using again the functional derivative approach, that our approximation scheme will be developed 
in section III. 



A. Definitions, equations of motion 

We work with creation-annihilation operators V-f, Vu for Wannier states of spin a =|, J. located at position ri, and, 
in the Heisenberg representation, imaginary time t\ . The space and imaginary time indices are abbreviated by arabic 
numerals. Furthermore, we use the Nambu representation, 

*t ( i) = (V{(i),^(i)) 

where the field operators obey the anticommutation relations 

{* a (l),*J(2)}j(ri-T 2 ) = J(l-2)(y ai/ 3 . (3) 

In this notation, the space part of Dirac delta functions are Kroenecker deltas: 5(1 — 2) = 6 TliT2 S (ri — T2) . When 
numerals are set in bold face, we mean just the space position (e.g. £(1,2)) and refer to the Schrodinger picture. 
Adding the convention that indices with an overbar are summed over space positions and, when appropriate, integrated 
over imaginary time from to /?, the Hubbard Hamiltonian takes the form 

H = -t (1,2) (4 (T) fa (2) + 4 (2) (1)) + (T) ^ (T) (T) Vi (T) 

with t (l,2) the hopping matrix elements. 

Since we will work in the grand-canonical ensemble, it is convenient to take H — fiN as the time evolution operator 
in the Heisenberg representation, with \i the chemical potential and N the number operator. The corresponding 
equations of motion for the field operators then are, 

S ( n - 75) W (2) = -Utf_ a (l) V-a (i) (i) (5) 



6 (n - 75) 4>\ (2) - u^_ a (l) (i) V£ (1) • (6) 
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(2) 



(4) 



d 

<9ti 



t(l,2) 



<9ti 



+ A £ ri , r¥ + *(l,2) 



We will also need the momentum space representation 



ik-ri t 



N 



e L c- 



k.tr 



and the pair operators 



with Fourier transform 



A(i) = v T (i)^(i) ; At (i)=^{ (i)^}(i) 



< c -k+ q a c k,T 



(7) 



(8) 



(9) 



These pair operators do not include the interaction potential in their definition. The pair operators obey the equations 
of motion, 



dA q = -1 \ - 

dr VJV^ 



£k + £-k+q - 2 Lu 



Ck,TC_k+q,| 



(10) 



5AJ, 



q _ 1 

, /at Z-^l 



£k + £-k+ q - 2 At 



'-k+q,| C k : T 



9r V^V k 

where the band dispersion e k is the Fourier transform of the hopping matrix elements 

-t(i,2) = ±-j2 elMri ~ r2)£ *- 



(11) 



(12) 



B. Green function and self-energy 



As in Ref. [3] , we work in the grand canonical ensemble in the presence of auxiliary source fields that are useful in 
intermediate steps of the calculations. The source fields are set to zero at the end. More specifically, we define the 
expectation value of a general time-ordered operator O by 



( T r [0]) @ = Z' 1 ({©}) Tr \ e 



with 

Z({&}) = Tr{e- f3 ("-^)T T [e 
where T T is the time-ordering operator while the matrix source field 

6(1,2): 

physically corresponds to Cooper pair sources, 



-*t(i)e(i, 2 )*(2) - 



* t (l)«(l, 2)*(2) 



0(1,2) 
6* (1,2) 



*t (l) (1, 2) * (2) = 9 (1, 2) </>{ (I) ^ (2) + 9* (1, 2) ^ (l) ^ (2) 
The Nambu Green's function is then a functional of the auxiliary field defined by 



G (1,2; {©}) = - <T r [*(!)*+ (2)]) e 



({T, 


"V T (1)V| (2)" 


[(T, 


~^(l)Vj (2) 



J / © 



<T r [V T (l)Vi (2)]) e 



^(l)V-l (2) 



(13) 



(14) 



(15) 



(16) 



(17) 
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The quantities that will appear at the end of the calculation will, in analogy with the Green function, be defined in 
zero source field, namely 



G (1, 2) = G (1, 2; 0) = - (T T [* (1) ¥ (2)] ) . 
Using the equations of motion (5) and following the algebra of Ref. [17], one finds that 

G- 1 (1, 2; {©}) = G - x (1, 2) - £ (1, 2; {©}) - (1, 2) 



(18) 



(19) 



where in Matsubara frequency (fc„ = (2n + 1) nT) and Fourier space (k), the non-interacting Green function takes its 
usual form, 



G (k,ik n ) 











(20) 



-ifc„-(£_ k -^i) J 

while the effects of interactions are contained in the self-energy matrix defined by 



S(1,3;{0})G(3,2;{©}) =U 



v>!(i + )v>i(i)^t a)^|(2) 

^(i+)V T (i)V>l(i)V4 (2) 



^(1 + )^(1)^t(1)^1 (2) 
^(1 + )Vt (1)^(1)^ (2) 



© 



(21) 



The notation 1 + indicates that the imaginary time is infinitcsimally larger than n, (or smaller for 1 ) . 



C. Self-energy and pair susceptibility 



The self-energy, as should be clear from the last equation, depends on four-point functions that may be calculated 
in different channels. For the repulsive Hubbard model, charge and spin fluctuation channels are dominant, so 
approximations for the four point functions are written down in these channels. However, in the nearest-neighbor 
attractive model away from half-filling, it is the pair fluctuations that are dominant, even in the paramagnetic state. 
As a preliminary remark, we can suggest how the self-energy will be related to the pair correlation function by making 
the key observation that the diagonal terms can be written as functional derivatives with respect to the 9 field. First, 
note that when 9 is set to zero in Eq.(21) it becomes diagonal and Sn simplifies to 



S n (1, 2) = -U (t t (1+) Vi (1) V>T (1) V>j (3)] ) G n X (3, 2) . 



(22) 



The operators located in the middle, tpi (1) ip-\ (1) , can be obtained from a functional derivative with respect to 
9* (1, 1) before this field is set to zero. Hence, taking the functional derivative and setting off-diagonal terms to zero 
afterwards, one is left with 



S 11 (l,2) = -C7 



<5G 21 (l+,3;{0}) 



59* (1,1) 



(23) 



©=o 



This is basically what we want. The self-energy will be written in terms of a response function in the particle-particle 
channel. 

To continue more generally, we step back and define the susceptibility matrix 



X(l,2,3,4;{©}) = 

where <5/<50 (3,4) is a matrix operator in Nambu space 

S ( 



6 



(50 (4,3) 



G(l,2;{©}) 



.50(3,4) 



56(3,4) 



SO* (3,4) 





(24) 



(25) 



such that 
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S® (3,4) 



0(1,2) = 5 (1-3) 5 (2 -4) I 



(26) 



where I is the identity matrix in Nambu space. Note that the susceptibility X is still a matrix in Nambu space 
with only two matrix indices (four matrix elements). With this notation, the two-point function that we will need to 
compute the self-energy from Eq.(21) is the following special case of Eq.(24): X (1, 2, 1 T , 1 T ; {©}) . Evaluating the 
functional differentiation explicitly, we have 



X (1,2,1* !*;{©}) 





se(i+,i+) 



50*(l-.l-) 







> T (l)Vj (2)" 


{-(T, 


(2)' 



e 



-<T T [^(1)^ (2)]> e 



— (T, 



(2) 




V'1( 1 )V'|(2)J) 35^7(^^1(1)^(2) 



(T T [^ T (l)Vi (2)]) 







© 



V>t(i)V>| (2) 

^(l)^(l-)^ T (l-)^| (2) 
^ T (l)^ T t (l+) ? A{(l+)^ T t (2) 
-F(l=F,l=F;{e})G(l,2;{©}) 
where we defined a function F that contains only the anomalous pieces of the full Green function 



V i t (l)Vi(l-)V'T(l-)V'i(2) 
^ T (l)^(l+)^(1+)^ (2)' 



F(l=F,l=F;{e}) 



— IT* 





',t „/,t 



Vl(i + )V T T (i + ) 



-(T T [V T (l-)Vi(l-)]) 




(27) 

(28) 

(29) 
(30) 

(31) 



The last term in the above expression for the pair susceptibility comes from functional differentiation of the denomi- 
nator in the definition of averages, as one can check with the definition Eq.(14) for the partition function, 





se(i+,i+) 







1 



Z({&}) 





Vj(i + )V>l(i + ) 



{T T [^(i-)^t(r)]) e \ -i 

J Z({&}) 



(32) 



Using rotational invariance and our general result for the self-energy Eq.(21), we finally find the following key 
relation between self-energy and susceptibility 



or 



S (1, 3; {©}) G (3, 2; {©}) = UX (l, 2, 1*, 1*; {©}) + UF (l* 1=F ; {©}) G (1,2; {©}) 
S (1, 2; {©}) - f/X (1,3, 1^, 1T ; {©}) G- 1 (3, 2; {©}) + UF (l=F, l=F ; {©}) 5(1-2)1 . 



(33) 
(34) 



D. Bethe-Salpeter equation for the three-point susceptibility and relationship between irreducible vertices 

and self-energy 

In this section, we derive the Bethe-Salpeter equation for the susceptibility using the functional derivative scheme 
[3] . This equation allows to define the two-particle irreducible vertex that plays for the susceptibility a role analogous 
to that of the self-energy for the Green function. In deriving the Bethe-Salpeter equation, we will recover the well 
known relation between self-energy and particle-particle irreducible vertices. Since the susceptibility is the functional 
derivative of G while the self-energy is trivially related to G _1 , it is natural to start from 

G(1,4;{0})G- 1 (4,5;{0}) =5(1-5)1 (35) 

and to take the functional derivative of this equation. One finds, 

[G (1,4; {©}) G 1 (4, 5; {©})] = (36) 



5 



(50(2,2) 



G(l,4;{©}) 



G- 1 (A,5;{@}) =-a+G(l,4 

-a_G(l,4) ( 



50* (2,2) 
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G- 1 (4,5;{©}) 



59 (2, 2) 



G- 1 (4,5;{0}) 



(37) 
(38) 



where we used the definitions 



0+ 











s 


1 


; = ( J 





Using the expression for the inverse Green function Eq.(19) and the chain rule to take into account the dependence 
of £ on {0} through [18] G, one obtains, 



SB* (2,2) 



G 1 (4, 5; {0}) = -5 (2 - 5) 6 (2 - 4) ct_ 



5£(4,5;{Q}) SG- e (6,7; {0}) 
<5G M (6,7;{0}) SB* (2,2) 



(39) 



SB (2,2) 



G- 1 (4, 5; {©}) = -5(2-5)5(2-4) a+ 



<5£(4,5;{0}) 5G^(6,7;{©}) 



*%(6,7;{e}) 



'(2,2) 



(40) 



The above equations simplify greatly when the functional derivative is evaluated in the normal, zero source field 
case, where particle number conservation implies the vanishing of anomalous correlation functions. We are left with, 



5G M (6,7;{&}) 



SB* (2,2) 



6G 2 i (6,7;{0» 



0=o -^Ai SB* (2,2) 



(41) 



©=o 



SG Te (6, 7; {0}) 



SB (2,2) 



<5Gi 2 (6,7;{0|) 
0=o ^ 2 50(2,2) 



(42) 



©=o 



Using the latter results in Eq.(37) multiplying it by G (5, 3) and then integrating over the point 5, we find, with the 
help of Eqs. (39) and (40), 



<50 (2,2) 



G(l,3;{0}) 



= cr+G (1, 2) a G (2, 3) + a G (1,2) er+G (2, 3) 



©=o 



5S (4,5;{0» 

+ " +G(M) 

+a " G(l ' 4) 5G 12 (6,7) 



5G 21 (6,7;{0» 



©=o 



59* (2,2) 
5G 12 (6,7;{0}) 



G(5,3) 



©=o 



©=o 



5B (2,2) 



G(5,3) 



(43) 
(44) 

(45) 



©=o 



Since we consider the normal phase, the off-diagonal Green functions vanish and we are left with only two equations 
that come from the diagonal components of the above matrix equation. The off-diagonal parts just tell us that two 
of the irreducible vertices vanish. Finally then, 



5G 21 (1,3; {©}) 



59* (2,2) 
<5G 12 (1,3;{©}) 



= G 22 (l,2 ) G n( 2,3 ) + G 22 ( 1 ,4)^|fll 
©=o 0G21 (0, 7) 



5G 21 (6,7;{0» 



©=o 



SB (2,2) 



©=o 



_ ftl(1 , J)ftlft , )+ai M^Wl 



59* (2,2) 
5G 12 (6,7;{0» 



Gn (5,3) (46) 



©=o 



©=o 



59 (2,2) 



G 22 5,3 



©=o 



(47) 
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E. An exact expression for the self-energy where low- and high-frequency behaviors are separated. 



The high-frequency limit of the self-energy is given by the Hartree-Fock result, as can be shown from sum rules. 
For latter purposes in our approximative scheme, it is useful to have at hand an exact expression for the self-energy 
where the high-frequency behavior appears explicitly. In this section, we derive such an expression in the case where 
the auxiliary field vanishes. 

First, let us recall the sum rules that fix the high-frequency behavior of the self-energy. In the absence of external 
field 



G(k 



,ik n ) — J 



cko A (k,w) 



2tt ik n — lu 

where the single-particle spectral weight A (k,w) for the Nambu Green function is given by, 



.(k,c)= J 



die 



iujt 



({ck, T (*).4,t} 







(*),<*,;} 
The high-frequency expansion of the Green function is then 



G (k A) = k I £ A ^ + jhl d ^ A M + • • 



(ik n ) 

where the frequency moments of A (k,w) are easily computed from equal-time commutators 



I 



2tt 



A(k,w) = 1 



Comparing with 



V -<n T ) 



G (k,ik n ) = (ik n I- (e k - A») a z - S (k,ifc n )) 1 

[(e k - jj) a z + S (k,ifc„)] + . . . 



(48) 



(49) 



(50) 



(51) 



(52) 



ik n [ik n ) 2 

it is clear, using spin-rotational invariance, that 



lim £ (k,ik n ) = U (n^) a z 



(53) 



(54) 



An expression where this asymptotic behavior is explicit may easily be obtained from the general formula for the 
self-energy specialized to zero external field, namely Eq.(34) with 9 = 0. Returning to our general discussion Eq.(23) 
will help understand the point we were making. In zero external field, we have 



S (1, 2) = UX (1,3, 1=F, IT) G 1 (3, 2) = -U -^^-^G (l, 3; {0}) 
or, looking only at the non-vanishing elements, 



e=o 



0^(3,2) 



En (1,2) =-U 
S 22 (l,2) =-U 



<5G 2 i(l+,3;{©}) 



60* (1,1) 
<5G 12 (1,3;{0}) 



Gii (3, 2) 



©=o 



50 (1+,1+) 



G^ 1 (3,2) 



(55) 

(56) 
(57) 



©=o 



Substituting the Bethe-Salpeter equations (46) (47) in the last two expressions we have the equivalent expressions, 
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E U (1, 2) = -UG 22 (1+ 1) S (1 - 2) - E/G 22 (1+3) ^' ( 2 _ { ® }) 



SG 21 (6,7;{©» 



©=o 



<S6>* (1,1) 



(58) 



©=o 



£ 22 (1, 2) = -UGu (1, 1+) «S (1+ - 2) - f/Gn (l, 4) 



<5£ 12 (4,2;{©}) 



<SGi2 (6, 7) 



SG 12 (6,7;{0» 



©=o 



66 (1+,1+) 



(59) 



©=o 



One may easily check that the terms proportional to 6 (1 — 2) in these two exact expressions are precisely the Hartree- 
Fock contribution, Eq.(54). 

The skeleton diagram representation of the Bethe-Salpeter equation (46) and of the self-energy Eq.(58) appear in 
Fig.l of the accompanying paper. This skeleton diagram representation is not necessary to understand the rest of 
the paper but may be useful for physical understanding. The non-perturbative approach developed in the following 
section does not directly correspond to the summation of an infinite subset of diagrams. 



III. A SYSTEMATIC NON-PERTURBATIVE APPROACH 



A brief reminder of how Hartree-Fock theory is derived in the functional derivative approach will motivate the two- 
particle self-consistent approach. At this first level of approximation (TPSC), one hasGW, in the presence of © 
and the corresponding irreducible vertices and susceptibilities are obtained from the functional derivative approach. 
The only unknown quantity, double occupancy (n^n^), may be obtained self-consistently by using what we will call 
the local-pair sum rule. The local-pair sum rule is a simple consequence of the fluctuation-dissipation theorem. An 
improved approximation for the self-energy, X^ 2 -*, will be found in section HID. We conclude with a discussion of an 
internal accuracy check that, as in the repulsive case, helps delineate the domain of validity of the approach. 



A. A perspective: Conserving approximations and Hartree-Fock theory 

In conserving approximations, the self-energy is obtained from a functional derivative of the Luttinger-Ward func- 
tional [19] that is computed from skeleton diagrams. Irreducible vertices for response functions are then obtained from 
appropriate functional derivatives. The Hartree-Fock approach is a special case of conserving approximation. The 
more standard way to derive the Hartree-Fock approach is to treat the general equation for the self-energy Eq.(21) 
in the presence of the auxiliary field as if Wick's theorem applied to the right-hand side of that equation. More 
specifically, Eq.(21) becomes, in the Hartree-Fock approximation 

E** (1,3; {©}) G HF (3, 2; {©}) = U ( ~f£ _f( F ft \} & ^ }) ) G HF (1, 2; {©}) (60) 

or 

^(Wfe)) = »(-||i}^'l -5$'e$<2>)) 4 < 1 - J > ' (61) 

The corresponding irreducible vertices for the Bethe-Salpeter equation governing the pair fluctuations, Eq.(43) are 
5Ef/(l,2;{©}) 52^(1, 2; {©}) 



«5Gf/(3,4;{©}) 



©=o «? 2 T(3,4;{0}) 



= I7J (1-2) £(1-3)5 (1-4) (62) 
©=o 



which leads to the simplest T-matrix type approximation. In the final calculations, we only need the self-energy in 
zero field. There it is purely diagonal and given by 



(1,2) = a z U( ni )S(l -2) . (63) 



Returning momentarily to the original equation that was approximated, Eq.(21), it is important to realize that the 
Hartree-Fock approximation satisfies an exact result. First note that the four point function becomes simple when 
point 2 becomes the same as 1. In that case, this quantity takes a simple form that, because of the anticommutation 
relations, depends on whether 2 — > 1+ or 2 — > 1~. More specifically, we have the exact results 
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S (l,3;(e } )G(3,lM6 ) )= I (l"'<";'l 1 »e (M^iiHBe) (M) 

E WW )GHl-:W).^ig<^. <„,(!) °,(D)e) ' (65 » 

In zero field, the difference of these two exact results is 

S(l,3)G(3,l+)-S(l,3)G(3,l-) = <r z U (n L ) . (66) 

Given what we have said in section HE about the exact high-frequency behavior of the self-energy, the reader will 
not be surprised to learn that the Hartree-Fock approximation does satisfy the relation Eq.(66). Indeed, this exact 
relation is sensitive only to the high-frequency limit of the self-energy, as may be proven by writing down explicitly 
the 1, 1 component of Eq.(66) in Fourier-Matsubara space as follows 

S ( k ' ifc ») U(ni) \( -ifc»o- c -»fc„o+\ (67) 

N Z^Z^\ ik n -(eu -u)-E(k Jk n ) ik n \ J K ' 



= U( ni ) (68) 



k ik n 



T 
N 



k ik n 



EE ^i>(e-«--e-'«*) 



In this expression, the first sum vanishes because we have added and subtracted a term that makes it convergent with- 
out the need for convergence factors e~ lfc ™° . Hence, only the last sum survives. The result is a direct manifestation 
of the anticommutation relations in the asymptotic behavior of the Green function. 



B. Two-particle Self-Consistency and irreducible vertex 

We will call the two exact results Eqs.(64) and (65) Tr[EG] for short cut. They are simply related to the potential 
energy (and hence to double occupancy), a crucial quantity for the Hubbard Hamiltonian. Furthermore, they can 
be considered as initial conditions for the true expression that defines the self-energy in the most general case, 
E (l, 3; {©}) G (3, 2; {©}) where point 2 has moved away from point one. Hence, in the first step of the approximation 
that we propose, we perform, as in conserving approximations, a Hartree-Fock like factorization in an external field, 
but we add the constraint that the factorization becomes exact when 2 — > 1 + or 2 — > 1~. More specifically, in analogy 
with the factorization Eq.(60) we start from the ansatz 

S«(1,3;{0})G«(3,2;{0}) =c^( ~ G S (1+ ' 1;{0}) G }( ! (1, 1; ' { ® }) | A ({©}) G^ (1, 2; {©}) (69) 
v ' ,x n v ' ,x n \ G$ (1, 1; {©}) -G${1, l+;{©})y 

where the matrix A ({©}) is obtained by requiring that the exact relations Eqs.(64) and (65) be satisfied [20]. Setting 
alternatively 2 — > 1+ or 2 — > 1~ in the last expression, and requiring equality with the respective exact Tr[EG] 
expression, we find, using spin rotational invariance, that in either case there is a unique solution for the off diagonal 
elements 

e8' ft * W) - v (Mi-» 1 »,ff(i.i;{e>)Mi-2) 

(»t>e (1 - »l)e - G « <L {e})Gg> (1. 1; {e}) 

with the analogous expression when the Nambu matrix indices 1 and 2 are inverted. It should be clear from this result 
that the functional dependence of A ({©}) on external field is only through G^ (1, l"* 1 ; {©}) and double occupancy 
( n l n l)@ ■ Nevertheless, this establishes a strong self-consistency relation between one and two-particle quantities that 
is absent from any standard diagrammatic approach. That is why we called this portion of the approach "Two 
Particle Self-Consistent" (TPSC) [11]. It is important also to note that the superscript (1) refers to what, in earlier 
publications, [12] wc called the zeroth step of the approximation. As we will see in a moment, at this level of 
approximation the diagonal part of E^ 1 ^ is a constant when = 0. Hence, in this limit, is equal to a bare 

propagator once the chemical potential is adjusted to obtain the proper filling. That is why we had referred to this 
level of approximation with superscript (0) in earlier publications. It is important to notice however that E' 1 ) (or 
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in the former notation) may have a dependence on external magnetic field, for example, that is absent in the 
non-interacting case. 

The particle-particle irreducible vertices appearing in the Bethe-Salpctcr Eqs.(47) and (46) are obtained from 
functional differentiation of the self-energy Eq.(70a) as in any conserving approximation. Since when = all off 

diagonal function such as G^ (1, 1) , G^ (1, 1) or 6 (n± (1 — n^)) @ /SG^ (1 5 1; {®}) an d the like vanish, we are left 
with 



6^(1,2;{®}) 



<5G«(3,4;{0}) 



5S«(1,2;{0}) 



©=o 



6G£>(3,4;{&}) 



©=o 



UppS (1-2) 6 (1-3)6(1 



U ^r^U d-2) 6(1 
(n L ) (l-n T ) 

4) • 



3)5(1-4) 



(71) 
(72) 



The irreducible vertex U pp is still local, as in the Hartree-Fock approximation, but, contrary to Hartree-Fock, the bare 
vertex is dressed. The function that dresses the vertex is simply related to double occupancy. It can be determined 
from the Bcthe-Salpeter equation itself by using the fluctuation-dissipation theorem, allowing us to close the system 
of equations. 



C. An approximate expression for E' 1 ' in the TPSC approach 

The self-energy entering G*- 1 ) in zero external field is diagonal. There is however an apparent paradox. Indeed, 
substituting the TPSC factorization Eq.(69) on the left-hand side of the two Tr[£G] equation Eqs.(64) and (65), all 
for = 0, seems to give two different solutions for the diagonal value of A ({0}) , and hence for the corresponding 
Let us use X^ 1 ) with a + or — index depending on whether they are the solution of Eqs.(64) or (65). In cither 

case, £W (1 - 2) is proportional to 6 (1 - 2) , so that using G$± (1, 1+) = (n T (1)) and G$ ± (1, 1") = -1 + (n T (1)) 

and the like for G^, one finds, 

E«+ (1, 2) = (1, 2) = (U- U pp (1 - n ; » 6(1-2) (73) 

(1, 2) = (1, 2) = U pp (n T ) 5(1-2) (74) 

From these results, keeping the filling fixed we have 

S« (1,3) G« (3,1+) - £ W (1,3) G W (3,1-) = a z U (n T (1)) (75) 

This suggests that the antisymmetric combination — is related to the high-frequency asymptotic behavior 
of the self-energy, as in the Hartree-Fock case. The symmetric combination on the other hand should not depend on 
convergence factors, as can be seen by using arguments analogous to those of Eq.(67). This quantity should then be 
a measure of the low-frequency behavior of the self-energy. With n = (nj.) + (n-f) the quantity 



\ ( E £> (1, 2) + (1, 2)) =*,(*- Ml^)) 6(1-2) 



(76) 



we expect, plays the role of a chemical potential shift, fj,^ — /io, with respect to the non-interacting value We use 
the notation fj,^ with a superscript (1) to note that this is the chemical potential corresponding to the self-energy 
We will give additional supporting evidence for this conjecture in section IV that discusses exact results satisfied by 
our approach and in Appendix A. In actual calculations the chemical potential used in the Green function occurring 
at this step of the calculation is determined from the condition G^ (1, 1 + ) = (nj (1)) with 

G$ (k, ik n ) -i (77) 

ik n - £k + M ( ' - s ii ( k , ik n ) 

(k, ik n ) in the combination fi^ — (k, ik n ) is a constant so that /j,^ — (k, ik n ) equals the non-interacting 
chemical potential yu - 

It is important to keep in mind the following. The approximate analytical expression Eq.(76) for and the 
corresponding chemical potential ^S 1 ^ are useful for analytical arguments and numerical estimates [21]. However, 
more fundamentally, the chemical potential is a thermodynamic quantity that can be computed in a consistent way, 
along with all other thermodynamic quantities. In diagrammatic methods, the Luttinger-Ward functional provides a 
systematic method to obtain such thermodynamically consistent results. [19] Since our approach is non-perturbative, 
another approach must be used. At this point, is a first estimate for the true chemical potential. A better 
estimate is jjS 2 \ corresponding to £( 2 \ the next level of approximation for the self-energy. 
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D. TPSC plus improved approximation for the self-energy 



Let us summarize what we have up to now. In this TPSC approach for the attractive Hubbard model, there is only 
one particle-particle irreducible vertex U pp Eq.(71). In addition, the two pair susceptibilities that we need, namely 



8G21 (2,2;{0» 



86* (1,1) 
5G 12 (2,2;{0}) 



= — (T, 



SO (1,1) 



= — (T, 



®=o 



^(1)Vt(1)^{ (2)V| (2) 
^{(1)^(1)Vt (2)V| (2) 



(78) 
(79) 



are quite easily related by the operation 2 <-> 1, and anticommutation ?/>j (2) ^ (2) = — i/>j (2) (2). So we define 



X P (1,2) 



SG21 (2,2;{0}) 



(1,1) 



SG 12 (1, 1; {©}) 



©=o 



'(2,2) 



(80) 



©=o 



The equality of these two functions is a reflection of the fact that for on-site pairing, the Pauli principle makes the 
triplet channel vanish. Hence, substituting the particle-particle irreducible vertex Eq.(71) in one of the Bethe-Salpeter 
equations Eq.(47) and expressing the result in terms of the ordinary Matsubara Green functions 



(81) 
(82) 



G; i} (1,2) =G«(1,2) 
G< 1J (1,2) =-G&(2,l) 

we find the pair susceptibility at the first level of approximation Xp^ (1,2), 

X« (1, 2) = G« (1, 2) (1, 2) - U„GW (1,4) X « (1, 2) G< 1} (1,3) . (83) 
In Fourier space and with the definitions, q = (q, iq n ) , g„ = 2nirT 1 k — (k, zfc„) , k n = (2n + 1) ttT and 

x^ ) (?)=§E G r ) ( fc +^ G i i) (- fc ) 



(84) 



for the particle-particle susceptibility that is irreducible with respect to U pp , the above equation will take the form 



4 1} (1) = 



1 + u„ x ^ (?) 



To close the set of equations, we need to solve for the irreducible vertex 

TJ _ rr ( n l C 1 ~ n T)) 
PP ~ <»i> (1 - »t> ' 

That may be done by using cither of the two exact limits 2^1+ and 2 — > 1 _ of Eqs.(78)(80) 

X « (l,l+)=<»l»T> 



X 



« (1,1") =l-n+<n i n T ) 



(85) 



(86) 



(87) 



(88) 



The left hand side of the latter two equations, that one could call local-pair sum rules, can be transformed into the 
following two sum rules when our approximation is used for the susceptibility 



X$ (1) 



4 1 )(l,l+) = (n i n T > = ^. 



-iq n 0~ 



(89) 



4 1} (M-)=l-«+<»l»T> = TfE 



Xlr (g) e - iqn Q+ 



N ql + U ppX t ] (?) 



(90) 
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We will discuss in section IV why, in our approach, both sum rules are consistent and so give exactly the same result 
for Upp. Monte Carlo simulations confirm that the pair susceptibilities calculated with Eqs.(86), (89), (84) and (77) 
are excellent approximations from weak to intermediate coupling. 

As we saw above, at this stage the self-energy X^ 1 ) entering the calculation is a constant determined in such 
a way that we have the proper filling. At high-frequency the exact limit of the self-energy is the Hartree-Fock 
result and the corresponding irreducible vertex is the bare interaction. In our section on exact results, we have 
found an expression for the self-energy, Eqs.(58)(59), that is the sum of two terms, a first one which is the high- 
frequency Hartree-Fock behavior and a second one that only involves the low frequency behavior of Green functions, 
of irreducible vertices and of susceptibilities. As we have seen above, our approximate expressions for these quantities 
are consistent and are low-frequency approximations. Hence, we substitute in Eqs.(58)(59) the result Eq.(83) for the 

susceptibility — 5G$ (2, 2; {&}) /59* (1, 1) = (1) 2) , and the corresponding results Eq.(71) for the irreducible 



vertex (1, 2; {©}) /5G$ (3, 4; {©}) 



©=o 

and Green function One is left with 



©=o 

£< 2) (1, 2) - Ud» (1, 1+) (5(1-2)- UUppG^ (2, 1) x « (1, 2) (91) 

Appendix B 1 discusses the relation to the corresponding formula in the repulsive case. Note that all the terms on 
the right-hand side of this equation, including the irreducible vertex U pp , are at the same level of approximation. In 
particular, the irreducible vertex U pp that appears explicitly, is the functional derivative of the field-dependent S*- 1 -* 

that enters and xp ■ This is crucial for the quality of the approximation [12,22]. This type of consistency is 
absent in many modern self-consistent treatments whose self-energy contains renormalizcd Green's functions but with 
only bare vertices. [5] Going to Fourier space and using the full expression for the susceptibility Eq.(85), we have 

i( 2 ) <U\ — TT^ TT—S^TT r<W I h l „\ 



S y> (k) = U ni ~U-^ Up P GY> (-k + g) % - - (92) 

iV q l + UppXlr(q) 



Spin rotational invariance gives us the result for down spins. Note that one of the vertices is bare while the other is 
dressed, contrary to the case where there is a Migdal theorem. 

The superscript (2) on the last expression for the self-energy indicates that it is the next level of approximation. To 

(9) 

improve the susceptibility calculation we would need the irreducible vertices corresponding to (1, 2; {©}) which 
we do not have. Hence the calculation stops at this level. Physically, the collective modes are less sensitive to details 
of the quasiparticles, so they can be computed first with simple Green functions. The self-energy on the other hand is 
sensitive to the collective modes (they have zero-frequency Matsubara contributions contrary to fermionic quantities) 
and hence we have to take these modes into account when we want a better approximation for the self-energy. 



E. Internal accuracy check 

Either by returning to the derivation Eq.(56) of xj 2 ^ (1, 2), which involves a four point function Eq.(78) related to 
the pair susceptibility Eqs. (83,80) and double-occupancy Eq.(87), or by starting from the Fourier space expression 

Eq.(92) and using the local pair sum rule Eq.(89) for xp 1 ^ an d the corresponding sum rule for Xil\ onc finds that 
Eq.(91) satisfies 

4 2) (1,2) (2, l + ) = §£ 4 2) (k) G{ 1} (k) e- ik "°- = U (n,n T ) . (93) 

k 

U (n^n-f) entering this equation is exactly the same as that computed from the local pair sum rule Eq.(89). This 
result is analogous to that found in the repulsive model. Note that G^ enters the above equation. An indication 

of the accuracy of our approximations may be obtained by checking by how much S| 2 ' (l, 2) G^ (2, 1 + ) differs from 
[/(n^n-f) obtained from the local pair sum rule. We have checked that there is at most a few percent discrepancy 
between both calculations, except in the pseudogap regime. Note that the chemical potential entering g| 2 ' must 
be obtained from the number conservation equation. We refer to Ref. [23] for a discussion of Luttingcr's theorem in 
this context. 

In any approximation for the many-body problem, full consistency for all correlation functions is unachievable. For 
example, in conserving approximations [4,5] one starts from diagrams for the Luttinger-Ward functional and for the 
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corresponding free energy. Then a self-energy and irreducible vertices are obtained from functional differentiation. 
These quantities may then be used in the Bethe-Salpeter equation to obtain the pair susceptibility. From that pair 
susceptibility, one can compute double occupancy through the exact result Eq.(87). The latter double occupancy is 
in general different from the one obtained from S (l, 2) G (2, 1+) since it does not contain the same set of diagrams. 
So if the double occupancy obtained from the susceptibility Eq. (87) is integrated over coupling constant to obtain the 
free energy, the result will in general be different from the original free energy. Conserving approximations are not 
self-consistent at the two- particle level. Other criticisms of these approaches appear in Refs. [12] and [22]. 



IV. EXACT RESULTS SATISFIED BY OUR NON-PERTURBATIVE APPROACH 



We briefly discuss exact relations and consistency requirements that are satisfied by our approach. Details of some 
of the proofs may be found in Appendix A. 



A. Sum rules on single-particle spectral weight 

Consider first the single-particle properties. These should be calculated with G*- 2 -* that contains the self-energy 
Eq.(92) and the corresponding chemical potential. One can extract the moments of the corresponding spectral weight 
A^ 2 ) (k,w) from the high-frequency expansion of G^ 2 ) (k, ik n ) in analogy with Eq.(50). In the self-energy S^ 2 ^ (k, ik n ), 
Eq.(92), the Hartree-Fock contribution appears explicitly so that the exact high-frequency limit Eq.(54) is satisfied. 
This means that the normalization and first moment of A^ 2 ) (k,w) satisfy the exact results Eq.(51) and Eq.(52). 



B. Sum rules on pair spectral weight 



Concerning two-particle properties, more specifically the pairing susceptibility, there are a number of exact results 
that our approach satisfies. First there are the two local-pair sum rules Eqs.(89)(90) that are a consequence of the 
fluctuation-dissipation theorem. The value of U pp obtained from either one of them is the same. This is demonstrated 
in Appendix A. 

The other exact properties we shall consider concern the pair spectral weight. Using the definition of the pair field 
A Eq.(8) and of the pair susceptibility Eqs.(80)(78)(79) the latter may be written as 

Xex(l,2) = (T T A(l)At(2)). (94) 

The subscript ex stresses that we are, for now, considering properties of the exact pair susceptibility. The Lehmann 
representation, and the periodicity in imaginary time, allow one to show that 

X«(q,^)=/^^ (95) 
J 7T u-iq n 

where the time Fourier transform of the pair spectral weight x"x (Qi u ) IS defined by 

X ' e ' x (q,i) = i([A q (i),At ( 0)]) . (96) 

From the latter definition and the equations of motion Eqs.(10)(ll) one can show that the quantity x"x (li^) obeys 
the following sum-rule 



/ 



^X'L (q, = < [A q (0) , A q (0)] ) = 1 - n (97) 



where n is the filling that is obtained from the single-particle Green functions entering the calculation of x"x (Qi^)- 
We show in Appendix A that our approximate expression for the susceptibility Eq.(85) satisfies this manifestation 
of the Pauli principle exactly, for all wave vectors q. That is why we can use either of the two local pair sum-rules 
Eqs.(89)(90) to find self-consistently the value of {n^n{) . 

Proceeding to the first moment of the pair spectral weight, it is shown in Appendix A, that 
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J ^-ux'L (q, w) = ^ E ( £ k + ^-k+q - 2 - |) ) (1 - (™ kT ) - (n-k+qi» (98) 
^E( £ k + e -k+q) (! - 2 Kt» - 2 - I) (1 - n) . (99) 



Like the previous sum rule, this result is valid for all wave vectors q. It is a generalization of the /—sum rule to 
the case of the attractive Hubbard model, a sort of off-diagonal /—sum rule. It is generalized Ward identity that 
relates two-particle quantities on the left-hand side with quantities obtained from the one-particle Green function on 



the right-hand side. At half-filling, /i — ^ = 0, where there is an exact canonical transformation from the attractive 
to the repulsive Hubbard model, the above result reduces precisely to the /—sum-rule for the repulsive case. Again, 
the above expression Eq.(99) relates a two-particle quantity, on the left-hand side, to a single particle property, 
on the right-hand side. Neither the pair susceptibility nor the single-particle Green function are known exactly in 
our approach. Nevertheless, as shown in Appendix A, when our approximation x^ f° r the pair susceptibility is 
substituted on the left-hand side of Eq.(99), and our expression for the corresponding single-particle Green function 
G^ 1 ) is substituted on the right-hand side, the equation is satisfied exactly as long as one uses 

M (1) =^o + |-^(l-n) (100) 

for the chemical potential appearing on the right-hand side. This is consistent with the fact that the quantities 
entering the right-hand side of this equation must pertain to the single-particle Green functions used in the calculation 
of x'p (q> w ) on the left-hand side. Since the chemical potential entering G^ should be fj,^ = /io+ X^ to obtain 
the correct filling, use of the approximation Eq.(76) for X^ leads to the above result Eq.(lOO). Recall however that 
Eq.(76) for X^ is not rigorous. Nevertheless, away from the renormalized classical regime, where the self-energy is 
weakly frequency dependent, the chemical potential obtained from this formula differs little from the one obtained 
at the second level of approximation X^ 2 ) or from Monte Carlo simulations in the weak to intermediate coupling 
regime [14]. The result Eq.(99) can also be considered as a type of consistency condition between one and two-particle 
quantities analogous to the relation between Tr (XG) and U (n^n-f) discussed in section HIE. 

C. Miscellaneous 

The Mermin- Wagner theorem is satisfied by our approach. That theorem states that classical fluctuation effects 
prohibit continuous symmetry breaking in two dimensions. The proof follows the steps of Appendix A. 3 Ref. [12]. 
The Kosterlitz-Thouless-Berezinskii [24] transition on the other hand involves algebraic order and large-scale vortex 
structures that are absent from the present approach. This transition is thus inaccessible to us. Far from the critical 
point on the other hand, one can show, following the steps Ref. [25] in analogy with the repulsive case, that there is 
Kanamori-Brucckner (quantum-fluctuation) type screening of U pp which is given by the approximate formula 

v,r ~ t^tu • A -£w<r^i: [*!?('>]'■ (101) 

Finally, it is important to notice that all calculations are done at constant density. In particular, the Green 
function G^ entering the calculation of X( 2 ) in Eq.(92) is evaluated at the same density as the final result G^ 2 \ This 
is motivated by the existence of Luttinger's theorem which states that the volume enclosed by the Fermi surface at 
T = depends only on density, not on interaction. It would be unphysical to iterate from X^ to X^ 2 ) starting from 
a G^ whose Fermi-surface associated singularities are at locations in the Brillouin zone that never intersect those 
of G^ 2 \ The constant-density constraint ensures maximum overlap. This point of view is motivated by Luttinger's 
approach. It is discussed further in Ref. [26]. Luttinger's theorem should be satisfied to a very good approximation 
in our approach, as for the repulsive model [12]. 

V. CONCLUSION 

In this paper we have presented a generalization of the approach developed in Ref. [12] to the attractive Hubbard 
model. We first established a number of exact results that form the basis of the approximation method that we 
introduced in section III. The first level of approximation (Two-Particle Self-Consistent) is based on a Hartree-Fock 
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like factorization ansatz for the self-energy in the presence of an external off-diagonal field, Eq.(69). That ansatz 
differs from the standard Hartree-Fock factorization Eq.(60) by the presence of a constant matrix A (0) that forces 
the factorization ansatz to reduce to the unfactorized four-point function in the special case where the later involves 
only density-density correlations (double occupancy) . That ansatz leads to the irreducible vertex for pair fluctuations 
simply through functional differentiation. The resulting irreducible vertex, given by Eq.(86), depends on double- 
occupancy, a quantity that may then be determined self-consistently using fluctuation-dissipation theorem derived 
sum rules Eqs.(89) and (90). Either one of these local-pair sum rules suffices to close the system of equations since 
they are equivalent. That exact equivalence is satisfied in our approach because the normalization sum-rule for the 
pair spectral weight, Eq.(97), is obeyed. This sum rule is a manifestation of the Pauli principle. 

The self-energy E^ entering the single-particle Green function at that first level of approximation is constant. The 
value of this constant is irrelevant for the calculation of the pair susceptibility since we work at constant filling, which 
means EW can be absorbed in the chemical potential. Nevertheless, we have argued that Eq.(76) should be a good 
approximation for the value of the constant self-energy at that first level of approximation since it follows from a 
consistency requirement between the ansatz Eq.(69) and the two possible values of Tr (EG) • In addition, Eq.(76) for 
E^ formally allows the first-moment (/ ^r^Xp (Qj 1 ^)) sum riue 011 the P arr spectral weight to be satisfied (Appendix 
A) . That sum rule is the off-diagonal generalization of the /— sum rule familiar from the particle- hole channel. 

The rough approximation that the self-energy is a constant suffices to obtain a good approximation for the low- 
frequency pair susceptibility since collective modes do not generally depend strongly on details of the damping of 
the underlying fermions. Details of single-particle damping do however depend strongly on the collective modes. It 
is possible then to improve our approximation for the self-energy by using our results for the collective modes in an 
exact formula Eq.(58) for £ where the high-frequency limit appears explicitcly. That gives us a better (cooperon- 
type) approximation for the self-energy Eq.(92). Clearly, that approximation does not assume that Migdal's theorem 
is satisfied since one of the vertices is dressed. When Migdal's theorem applies, both vertices are bare. As in previous 
work, [12], one can use the difference between X^ 2 ' (l, 2) G^ (2, 1+) and U (n^n-f) as a check on the level of accuracy 
of approach, as discussed in Sec. HIE. Formal comparisons with other approaches are presented in Appendix B. 
Note that our approach is in the SO (N — > oo) universality class, by arguments similar to those that apply to the 
repulsive model. [27] Hence, the Mermin- Wagner theorem is satisfied as it should, but algebraic long-range order of 
the Kosterlitz-Thouless type is beyond the accuracy of any microscopic theory that does not use renormalization 
group arguments to reach the long wavelength limit and treat the SO (2) symmetry exactly. 

In summary then, in its simplest form, our generalization of Ref. [12] to the attractive Hubbard model is expressed by 
the three simple equations, Eqs.(86),(89) and (92) plus the constant density contraint that determines the chemical 
potential appropriate to the level of approximation. Extensions of this approach have also been proposed. [28] 
Questions related to thermodynamic consistency and calculation of other response functions will be presented in a 
later publication. In the accompanying paper, our approach is compared in detail with results of Quantum Monte Carlo 
simulations. In addition to achieving quantitative agreement with simulations, this approach predicts the appearance 
of a pseudogap in the single-particle spectral weight when the pair fluctuations enter the renormalized classical regime. 
This is qualitatively different from the results obtained from self-consistent T — matrix approximation, or FLEX- type 
approaches. The role of the low space dimension, the pair correlation length and single-particle thermal dc Broglie 
wave length, and more generally the mechanism for the opening of this pseudogap, have been discussed in detail in 
Ref. [12]. 
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APPENDIX A: SUM RULES AND HIGH-FREQUENCY EXPANSION FOR THE PAIR 

SUSCEPTIBILITY 

In this appendix we give more details on the derivation of the results presented in section IV. Let us begin with the 
high-frequency expansion of the pair susceptibility Xp (q, iQn) for any approximation that has a spectral representation 
such as Eq.(95). From that spectral representation, one obtains, 
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/duj ( 1 \ f duj ( 1 \ 



(Al) 



Let us now consider the exact susceptibility \ex (q, *<7n)- The moments of \" x (q, u) that appear as coefficients of the 
expansion in powers of (iq n ) may be obtained as follows. The first one follows for all values of q from the equal-time 
commutator in Eq.(97). The next coefficient, / ^wXex (q, ^) , follows from 



/ ^V4(q, W ) =<[!/ X L(q,-) 



t=o 



- [A q (t) , At (0)] 



t=o 



dr 



A q (r),At(0) 



(A2) 
(A3) 
(A4) 



T=0 



The latter equal-time commutator may be computed after ^A q (r) is rewritten with the help of the equation of 
motion, Eq.(10), leading to the result Eq.(99). 

Since our approximate expression Eq.(85) for the pair susceptibility admits a spectral representation, its moments 
may be obtained from the high-frequency expansion in powers of (iq n ) , by analogy with the method in Ref. [12]. 

From our approximate formula Eq.(85) xi (?) 1 
that 



X^ir (?) 1 + U pp and the large iq n expansion of x\r (?) one nn( is 



(i) 



lim iq n Xp 1] (q, iq n ) = . lim HnXt-' i( ln) = - (1 - n) 



(1) 



(A5) 



in agreement with the exact result Eq.(97). Note that the large iq n limit of x$ (q> *<Zn) must be taken after the sum 
over fcrmionic frequencies in Eq.(84). 

The first moment (off-diagonal / sum-rule) is given by 



- / — Mv"' 1 ' 



lo X ; W (q,w) = . lim 

7T i<j„^oo 



= lim 



(*9«) 2 X p 1} (q, iq n ) + iq n (1 - n) 

fan) 2 Xi^ (q, ^n) (l - J/ppXir 5 (Qj *9n)) + (1 - Tl 

Substituting the exact results for the high-frequency expansion of x\l (q> i one obtains, 

| (q, w) = ^ ]T (e k + e- k+q ) (1 - 2 (n kT » - 2/i (1 - n) + <7 PP (1 - n) 1 



(A6) 



(A7) 



In this expression, (rikf) is computed from G^, hence is a Fermi function, and fj,o = /i^ + X^ 1 ' enters the Green 

function G^ from which Xil^ is computed. Since the self-energy EW entering G^ is a constant, //o coincides with 
the non-interacting chemical potential appropriate for the filling we are considering. Comparing with the exact 



(i) 



Mo 



result Eq.(99) we see that the chemical potential at that level of approximation should be given by fj, 
(U — Upp (1 — n)) /2 which coincides with our proposed approximation Eq.(76). Numerically, we have checked that 
this approximate chemical potential is quite close to the chemical potential // 2 ) obtained with E' 2 ) and that the latter 
in turn is close to those obtained from Monte Carlo simulations [14], as long as we are far from the renormalizcd 
classical regime where S^ 2 - 1 acquires a strong frequency dependence. 



APPENDIX B: COMPARISON WITH OTHER APPROACHES 



We first comment on the connection to the formalism for the repulsive model. This leads us then to a discussion of 
the widely used Self-Consistent T— matrix approximation and of another "mixed" approach that has been extensively 
applied recently. [29] 
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1. Mapping to the repulsive model and comparisons with Self-Consistent T —matrix approximation 



At half-filling, the Lieb-Mattis canonical transformation 




(Bl) 
(B2) 



with Q = (7r, 7r) maps the repulsive into the attractive Hubbard model. The same canonical transformation maps 
spins S, density p and pairing operators into each other as follows 



For a chemical potential different from U/2 (half-filling), the repulsive model maps into the attractive model at half- 
filling but in a finite Zceman-coupled magnetic field. The approach presented here for the attractive model would, 
at half-filling, translate into the transverse-channel calculation for the repulsive model [22]. In that U > case, Ref. 
[30] presents only the longitudinal channel calculation. The analog calculation for the attractive Hubbard model 
would have lead us to two irreducible vertices. One vertex would have been for the charge fluctuations. These are 
related to pair fluctuations by the SO (3) symmetry of the model a n = 1, which implies that the corresponding 
irreducible vertex there is also U pp . The other vertex would have been for the non-singular spin-fluctuation channel. 
The two vertices would have appeared in a self-energy formula that would replace Eq.(91). However, at n = 1, the 
best self-energy formula for the attractive Hubbard model would be obtained from the canonical transformation of 
that presented in Ref. [22] which preserves crossing symmetry. For problems sufficiently far away from half-filling 
(Tx <C /x) however, the pair fluctuations suffice. [31] 

In the repulsive model case, we have presented general analytical arguments [12] as well as detailed comparisons 
between: Monte Carlo simulations [12,22], our approach, and self-consistent Eliashberg type approaches (such as 
the Fluctuation- Exchange approximation) . Most of our general criticism concerning self-consistent approaches in the 
repulsive case apply to self-consistent T — matrix plus fluctuation exchange approaches in the attractive case. 

More specifically, one of the key qualitative differences between our approach and self-consistent approaches is that 
the latter do not predict the existence of a fluctuation-induced pseudogap in the single-particle spectral weight in two 
dimensions. We have demonstrated at length, through comparisons with Monte Carlo simulations [22,32] and with 
physical arguments [12], that this is incorrect in both the repulsive and the attractive cases. 



An alternate approach based on computing the irreducible susceptibility with one bare and one dressed Green 
function has been extensively used lately [29] . More specifically, in this approach 



S z (Q + q) 
S + (Q + q) 
S- (Q + q) 



p(Q + q) 

-At (q) 
-A(q) 



(B3) 



2. The GG (0) approach 



Xp = 



(B4) 



where 




(B5) 



k 



with the self-energy entering G given by 




(B6) 



(B7) 



Since there is a single chemical potential for the two Green functions and G®, there are two different expressions 
for the occupation numbers operators (n^ and n£) and for the corresponding fillings, (n and n°). 
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A positive aspect of this approach is that it exhibits consistency between one and two particle properties in the 
sense that the exact result 

^ £ S A (k) G, (k) = U (mn T ) (B8) 

k 

that follows from the equation of motion, is satisfied exactly by the above approximate scheme. Indeed, starting from 
the approximate formula Eq.(B6) for the self-energy, and using Eq.(B5) for the susceptibility, one finds that 

J Y S A (k) G i (k) e- ife "°" = U J V - X r p r L g) , r e-^"°" (B9) 
AT A ^ 1 + C/xp (?) 

= £7 (A* A) = t/(n T «l> (BIO) 

One goes from the second to the last line using the fluctuation dissipation-theorem. Hence, in this approach, the 
double occupancy obtained from single particle quantities (namely £j.Gj,) is exactly the same as that found from the 
pair susceptibility, a two-particle quantity. 

On the negative side, the spectral weight corresponding to the susceptibility Eq.(B4) does not satisfy the sum 
rules on the first two moments discussed in Eqs.(99) and (97). To show this, we first need a few sum rules on the 
single-particle spectral weight. Following the steps in Appendix A of Ref. [12], the high-frequency expansion of G, 
with S given by Eq.(B7), gives the following sum rules for the corresponding spectral weight A (k,w) 

^A(k,o/) = l (Bll) 



/ 



puA(\i,uj)=£ k -fi + Un _ a . (B12) 



Also, the equation of motion for G gives 

duj 1 

—uf (w) A (k, u>) = - ( £ k - M) «k, CT + U (n T n ; ) . (B13) 



/ 



2tt " v ' x ' ' N 

k 



The above sum rules are valid for both G and G^-°\ In the latter case however, we take the Fermi function for the 
occupation number and U = on the right-hand side of the above equations. 

We are now ready to check the sum rules for x p . Using the spectral representation for G in the expression for the 
susceptibility Eq.(B5), 

and performing the high-frequency expansion after the summation over ik n , one obtains 

lim {iq n Xp(ci,iq n )) = -1 + rn+n* (B15) 

with the help of the normalization Eq.(Bll) and first moment sum rule Eq.(B12). This should be compared with the 
exact result (—1 + n) found in Eq.(A5). The difference between the two fillings and n° is a measure of how much 
this approximation violates the Pauli principle. 

Pursuing the large iq n expansion and using spin rotational invariance and the sum rules Eqs.(Bll) to (B13) on the 
single-particle spectral weight, one obtains for the first moment of the pair spectral weight (off-diagonal /—sum rule) 



/ 



duo 
— UX' V 



=^E ( £ k + e-k +q -2^-^jj (l-(n kii )-(n°_ k+q , T )) 

-Un l -U((n m )-n l n l ) + 2Un ] n l ] . (B16) 
Even with nt = n k there are deviations from the exact result that are linear in U. 
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